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Density matrix perturbation theory [Phys. Rev. Lett. 92, 193001 (2004)] provides an efficient 
framework for the linear scaling computation of response properties [Phys. Rev. Lett. 92, 193002 
(2004)]. In this article, we generalize density matrix perturbation theory to include properties 
computed with a perturbation dependent non- orthogonal basis. Such properties include analytic 
derivatives of the energy with respect to nuclear displacement, as well as magnetic response com- 
puted with a field dependent basis. The non-orthogonal density matrix perturbation theory is 
developed in the context of recursive purification methods, which are briefly reviewed. 
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I. INTRODUCTION 

During the last decade a new computational paradigm 
has evolved in electronic structure theory, where no crit- 
ical part of a calculation is allowed to increase in com- 

Linear scaling electronic structure theory extends tight- 
binding, Hartree-Fock, and Kohn-Sham schemes to the 
study of large complex systems beyond the reach of con- 
ventional methods. In general, conventional methods 
have three computational bottlenecks: (1) construction 
of the Hartree-Fock, tight-binding, or Kohn-Sham Hamil- 
tonian, (2) solution of the self-consistent-field equations 
to obtain the ground state, and quite often, (3) the eval- 
uation of response properties. Here we will focus on the 
last problem, including calculations of both linear and 
higher order non-linear response. The main purpose is 
to present a non-orthogonal generalization of the den- 
sity matrix pert urbation theory, recently introduced by 
the authors 2^, which was used for linear scaling com- 
putation of static electric polarizabilities by perturbed 
projection p^ . 

A non-orthogonal generalization of the density matrix 
perturbation theory is important because it includes ba- 
sis set dependent perturbations in the overlap matrix, 
i.e. the basis function inner products, when using non- 
orthogonal local orbitals. In this case, the perturbed 
Hartree-Fock or Kohn-Sham eigenvalue problem may be 
expressed in the generalized form 

{R + R')^,^e,{S + S')cl),, (1) 

with both the effective Hamiltonian H and the overlap 
matrix S modified by the perturbations H' and S', respec- 
tively. Such perturbations are encountered when using a 
basis of atom centered Atomic Orbitals and computing 
geometric energy derivatives js^ , and also when cal- 
culating magnetic response properties with field depen- 
dent Gauge Including Atomic Orbitals. 

Previously, non-orthogonal generalizations of linear 
scaling methods for calculation of the density matrix 



have been developed by including the overlap matrix 
as a metric tensor in operator products [23, 113, 13 ■ 
Non-orthogonal density matrix schemes for solving the 
coupled-perturbed self-consistent-field equations [33, 
I41I ] have also been put forward, which pose density ma- 
trix derivatives implicitly through a set of commuting 
Sylvester-like equations 42]. In this article, we present a 
non-orthogonal generalization of density matrix pertur- 
bation theory, based on an explicit recursive expansion 
of the non-orthogonal density matrix and its derivatives. 
This theory provides a framework for extending linear 
scaling response calculations to properties with a non- 
orthogonal basis set dependence on the perturbation. 

The article is outlined as follows: First we give a brief 
review of the orthogonal formulation of density matrix 
purification and density matrix perturbation theory. Op- 
erators represented in an orthogonal basis set are de- 
scribed by italics (P). Thereafter we generalize the de- 
scription to non-orthogonal representations, where the 
matrices are described by normal letters (P). The central 
result is the non-orthogonal density matrix perturbation 
theory in section Fill Bl A simple example is given in de- 
tail for the expansion of the interatomic pair interaction 
of the diatomic HJ molecule up to fourth order. 



II. ORTHOGONAL DENSITY MATRIX 
PURIFICATION AND DENSITY MATRIX 
PERTURBATION THEORY 

A. Orthogonal purification 

Linear scaling electronic structure theory is based on 
the quantum loc ality (or nearsightedness) of non-metallic 
systems [13, [H, [3[4l, [i^. In a local basis, this local- 
ity is manifested in an approximate exponential decay 
of density matrix elements with inter-atomic separation. 
Under these conditions, and using a sparse linear algebra 
with the dropping of numerically small elements below a 
threshold t, the number of non-zero entries in the density 
matrix scales asymptotically linearly, 0{N), with system 
size N. This sparsity is used to achieve an 0{N) com- 



plexity in an iterative construction of the density ma- 
trix by operating only with sparse intermediate matri- 
ces. Several of these techniques are based on the Fermi 
operator relation between the density matrix P and the 
Hamiltonian H taken at T = 0, 



(2) 



given by the step function 9 (spectral projector), with 
the step formed at the chemical potential /i. The chemi- 
cal potential determines the occupied states via Aufbau 
filling. 

One approach to constructing P is thro ugh exp ansion 
of 9 using the Chebychev polynomials Tk |47Ll48ll49t : 



(3) 
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FIG. 1: The iterative expansion of the step function using 
second order trace correcting purification, Eq. 



With the two-term recurrence 



Tk+i{X) = 2XTk{X) - Tk-i{X), (4) 

the computational cost for evaluation of Eq. scales 
as 0{p) with polynomial order p of the Chebychev ex- 
pansion. However, this cost can be reduced to 0{y/p) 
by using a hierarchical summation of polynomial terms, 
rather than the two-term recurrence relation [iM. l50l| . 
While this reduces the computational cost, it demands 
more intermediate memory to store temporary matrices. 
In either case, the order of the polynomial approxima- 
tion in Eq. (|3J), P, must be kept fairly low, leading to 
problems with incompleteness. An incomplete Cheby- 
chev series results in Gibbs oscillations, which are high 
frequency ripples that form about the step. These oscilla- 
tions can be reduced by applying Gibbs damping factors 
to the Chebychev polynomials as in the Kernel Polyno- 
mial Method However, this reduces the slope 
of the step function and a higher order expansion is nec- 
essary to resolve the step at the chemical potential [49l |. 

Alternatively, density matrix purification is a recursive 
app roach to spectral projection [23, IH, |^ [s^, IH, 
I54L I55I Is^ that approximates the matrix step function 
projection as 

P = 9{fil -H)^ lim F^{Fn-i{...Fo{H)...)). (5) 

n — >oo 

Initiating a purification sequence is the linear transform 
Fq{H), which normalizes the spectra of H to [0, 1] in re- 
verse order. The functions F„ {n > 0) are typically low 
order, monotonically increasing polynomials, with fixed 
points at and 1. Each purification polynomial Fn grad- 
ually shifts the eigenvalues of the approximate interme- 
diate density matrix X„ to for unoccupied states and 
to 1 for occupied states as 

X„+i = F„(X„) = . . . = F„{...FoiH)...) , (6) 

creating a successively more "purified" intermediate 

Xn+l- 



Purification has several advantages. First, a p order 
truncated purification sequence can be developed with 
a complexity of O(logp). For example, with quadratic 
polynomials, a, p ^ 10^ order expansion can be reached 
in only 30 iterations. Also, because the polynomials Fn 
are monotonically increasing, so to are the corresponding 
purification sequences, regardless of the degree of incom- 
pleteness. The Gibbs oscillations resulting from trun- 
cation of the Chebychev series are therefore avoided and 
the application of damping factors is no longer necessary. 
Figure n illustrates the typical behavior using second or- 
der trace correcting purification as described below. 

To achieve a linear scaling with purification, threshold- 
ing is applied after each recursive expansion step to the 
intermediates X„, removing elements below a tolerance 
r 10""* — 10~^). This often leads to a substantial in- 
crease in computational efficiency, but also to an accumu- 
lation of numerical error with each recursive purification 
step. This error involves corruption of the eigenbasis, 
which at first increases exponentially ^^51] . However, this 
error accumulation disappears as the eigenvalues of X„ 
approach or 1. Since the number of purification steps 
necessary to reach convergence scales with the logarithm 
of the inverse band gap, the method is stable and the to- 
tal accumulated error is well controlled. At convergence 
the density matrix error scales linearly with the thresh- 
old T and the error in total energy decreases quadratically 
with decreasing r [sil . 

Density matrix purification methods differ in the way 
the purification polynomials Fn{Xn) are chosen. In grand 
canonical schemes |23,|2a,|53 the initial linear normaliza- 
tion Xi — Fq{H) shifts the eigenvalues such that all oc- 
cupied eigenvalues are in [c, 1] and all unoccupied eigen- 
values in [0, c], where c is some predefined number (typi- 
cally c = 0.5). Thereafter a fixed purification polynomial 
with infiection point at c is used, which shifts eigenval- 
ues above (below) c to 1 (0) . At convergence the correct 
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occupation is therefore reached, with 

Tr(P) = Km Tt{X^) = N,. (7) 

n — 'oo 

The problem with this approach is that it requires prior 
knowledge of the chemical potential fi, which has to be 
shifted to the inflection point c in the initial normal- 
ization Xi — Fo{H). To avoid the problem with an 
unknown chemical potential Falser and Manolopoulos 
(PM) devised a canonical purification scheme [2^ with 
the purification polynomials chosen such that the trace, 
i.e. the occupation, is preserved in each purification 
step. By choosing the initial normalization such that 
Tr(Xi) = Tr[Fo(i7)] = the PM scheme automatically 
converges to the correct density matrix, without prior 
knowledge of the chemical potential. The problem with 
this method is that it has a very slow convergence at high 
or low occupation [20ll2^ . A solution to this problem was 
give n by the introduction of trace correcting purification 
l2al described below. 



1. Second order trace correcting purification 

Trace correcting purification [26l Isil l55l | is an efficient 
approach to density matrix purification at both high and 
low occupation and does not require knowledge of the 
chemical potential. In trace correcting purification the 
polynomials F„(X„) correct the trace Tr(X„) and ex- 
pand the step function simultaneously. At convergence, 
the correct occupation of the density matrix is reached, 
such that Tr(F) = lim„^oo Tr(X„) = Ne. The simplest 
and most memory efficient form is the second order trace 
correcting algorithm |26| . given by 

S^max ^ min 



Xn+1 — Fn{Xn) — Xn + (Jn{I — Xn)Xn, (9) 

where 

an = sign(Arc - iVocc), (10) 



iVocc = Tr(X„), (11) 

and 

P = lim X„. (12) 

n — 'oo 

The constants £max and £min are upper and lower esti- 
mates of the spectral bounds of H, given for example by 
Gersgorin bounds 0|. The initial normalization Fo{H) 
thus transforms all eigenvalues of H to the interval [0, 1] 
in reverse order. The function sign(x) denotes the sign 
of a;. It is -1-1 if X > 0, otherwise it is —1. The trace 



correcting recursion in Eq. is equivalent to previously 
published versions |2^ |2^ i2^J , but is presented here in a 
form more closely related to an efficient implementation. 

In the case of aggressive thresholding and high or 
low occupation, the eigenvalues of X„ can sometimes 
be pushed out of the domain of guaranteed convergence, 
[0, 1]. To enhance stability under loose thresholding, we 
alternate the sign of the trace correction in every step, 
with (7„ = — cTn-i as convergence is approached, typically 
when Tr[(/ - X„)X„] < 0.1. 

B. Orthogonal Perturbation Theory 

The main obstacle in formulating a density matrix per- 
turbation theory based directly on the relation between 
the Hamiltonian and the spectral projector, given by 
Eq. |(2J), is the discontinuous, non-analytic nature of the 
step function. Difficulties with this discontinuity are fur- 
ther amplified when considering direct expansion of the 
projector. At finite temperatures close to zero, we could 
use the analytic Fermi-Dirac function, but this involves 
the computation of matrix exponentials and requires the 
chemical potential a priori to high precision. However, 
purification methods furnish a recursive, analytic, mono- 
tonically increasing and highly accurate representation 
of the step function that does not require prior knowl- 
edge of the chemical potential. The fundamental idea 
behind our approach is that this representation can be 
used in a direct variation of the density matrix with re- 
spect to a perturbed Hamiltonian, where perturbations 
in H can be carried through at each level of purifica- 
tion, either exactly or to finite order (2^. At finite order, 
this theory provides a framework for the computation of 
density matrix derivatives in the A'^-scaling computation 
of adiabatic response properties by perturbed projection 
|29l l57| . At infinite order, i.e. with an exact expansion, 
the method can be used for efficient quantum embedding 
of local perturbations . 

1. Exact expansion (infinite order) 

Assume a perturbation in the Hamiltonian, 

H^H'-^^+H'. (13) 
The recursive expansion of the density matrix 

F„(F„_i(. . . F,{Fo{H^"^ + H')) . . .)) (14) 
generates the corresponding perturbed sequence, 

Xn — Xn + A„, {^5) 
Xn+1 = Fn{Xn), 

where X^^ is the unperturbed sequence generated from 
X^^l^ = Fnixi°^) with = H^°\ and the initial per- 
turbation Aq = H'. The perturbed orthogonal density 
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matrix is given by 



P = P(°^ + lim A, 



(16) 



with the purification differences 

A„+i = F„(X„ + A„) - F„(X„). (17) 

Combined with the second order trace correcting purifi- 
cation, Fn in Eq. we have the foliowing recursive 
scheme for infinite order, orthogonal density matrix per- 
turbation: 



Ai H /(^max ^min); 



2A„ - A„} - A2 , iVocc < iVe, 



(18) 



(19) 



where we use the anti-commutator notation {A, _B} ~ 

AB + BA and the occupation iVocc = Tr(Xn''''). Since 
the differences A„ change quadratically in each iteration, 
the expansion order is in practice infinite at convergence 
and therefore exact. For insulators the computational 
cost scales linearly with the size of the perturbed region 
0(7Vpert.) since the recursion only involves terms with 
the response factors A„. For a local perturbation the 
computational cost is therefore independent of system 
size, i.e. it scales as 0(1) p^ . 

The perturbation theory is grand canonical since the 
expansion of the perturbation is performed at a fixed 
chemical potential determined by the unperturbed (or 
the perturbed) system. For sufficiently large perturba- 
tions, states may cross the chemical potentials /z. In this 
case lim„_+oo Tr(A„) ^ and the system is no longer 
neutral. 



The occupation iVocc = TT:[Xn^). The density matrix 
derivatives are given by 



1 



i! 5A'= 



= P(™) = lim Xl' 



A=0 



such that 



P = P(") + AP(1) + A2p(2) 



(23) 



(24) 



These equations provide an explicit and rapidly conver- 
gent algorithm for the computation of the density matrix 
response to high order in the expansion parameter |28l| . 
In addition, the formalism presented here is remarkably 
simple and can be easily extended to multiple indepen- 
dent perturbations (29l.l57i|. 



III. NON-ORTHOGONAL DENSITY MATRIX 
PURIFICATION AND DENSITY MATRIX 
PERTURBATION THEORY 

In a non-orthogonal representation, the Hartree-Fock 
or Kohn-Sham equations may be posed as the generalized 
matrix eigenvalue problem. 



(25) 



where the overlap matrix S is a matrix of basis function 
inner products. 

In the following, a non-orthogonal density matrix pu- 
rification algorithm is developed. Then, a non-orthogonal 
density matrix perturbation theory is introduced that ad- 
mits simultaneous perturbations in both the Hamiltonian 
and the overlap matrix. Normal letters (H) are used 
to distinguish the non-orthogonal representation from 
the corresponding orthogonal representation denoted by 
symbols in italics {H). 



2. Finite perturbation expansion 
Assume a perturbation expansion of the Hamiltonian, 

H = + Ai/(i) + A2i/(2) + . . _ (20) 

This perturbation generates the corresponding perturbed 
sequence 

X„ = + AX(i) + \^X^) + . . . , (21) 

where the separate m"^ order perturbations X^™^ can be 
collected order by order. Using the second order trace 
correcting scheme and keeping terms through order m 
in A at each iteration, the following explicit recursive se- 
quence is obtained for m — TOmax, "^max — 1, . . . , 1, 0: 

{m 
^XW^i™-^), iVocc>^e 
(22) 
i=0 



A. Non-Orthogonal Purification 

In the non-orthogonal case, the necessary criteria de- 
termining the density matrix are 

STH - HPS = 0, 

Tr(PS) = TVc, (26) 
P = PSP, 

together with Aufbau filling, i.e. occupying the iVc lowest 
eigenstates. 

Following normalization, Xi — Fo(H), purification pro- 
ceeds as in the orthogonal case, but with the minor ad- 
dition of the metric S to each (previously orthogonal) 
matrix-matrix multiplication, i.e. 

Z = XY Z = XSY. (27) 

For the second order trace correction purification, Eq. 

P„(X) =X2 ^ F„(X)=XSX, 

Fn{X) =2X-X2 -> F„(X) = 2X-XSX. 
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This preserves the covariant (or contravariant) form after 
each purification step [ssj . 

The most challenging aspect of non-orthogonal purifi- 
cation is obtaining an initial normalization Xi, which 
must obey the commutation relation 

SXiH-HXiS = 0. (29) 

With this normalization, commutation is automatically 
preserved as long as 

SF„(X„)H - HF„(X„)S = SX„H - HX„S. (30) 

This is true for all non-orthogonal purification polynomi- 
als, which have the form 

t 

F(X) = ^afcX(SX)^ (31) 

k=Q 

where the are polynomial expansion coefficients. 

There are a number of options for initiating the non- 
orthogonal purification, three of which are 

Xi=Fo(H) = a(/3S-i-S-iHS~i), (32) 
Xi ==Fo(H) =a(H-i +/3S~i), (33) 
Xi =Fo(H) =a(H-/3S)-\ (34) 

where a and f3 are chosen to map the eigenvalues into 
[0,1] in reverse order. Ideally, an efficient choice of nor- 
malization concentrates all unoccupied states near 0, and 
all occupied states near 1. 

The first choice, Eq. H32|l . is analogous to the initial 
guess suggested by Falser and Manolopoulos |23| for their 
non-orthogonal grand canonical purification scheme, and 
amounts to a shift and linear rescaling of eigenvalues. 
This linear rescaling can lead to a small renormalized gap 
in [0, 1] for the case of low occupation in the large basis 
set limit. The second case, Eq. (|33|) . is of minor interest 
since it involves calculations of both and S~^. The 
third normalization, Eq. H34|) , is the most interesting and 
useful of the initializations. With a — 1 and (3 — emin^ 1, 
where Emin is a lower bound of the eigenvalues Si in Eq. 
(|25|l . this is a Green's function, 

G(/?) = (H-/3S)-\ (35) 

which provides the correct normalization. In this nor- 
malization, the unoccupied states are mapped to as 
l/si. For large basis sets, with a low fractional occupa- 
tion, this amounts to a rescaled band gap on the interval 
[0, 1] that is larger relative to the band gap given by the 
linear rescaling. Since the number of iterations needed to 
reach convergence scales with the logarithm of the inverse 
band gap [26] the Green's function initialization can be 
expected to be more efficient in the large basis set limit. 

1. Computation and refinement o/Xi = G(/3) 

The Green's function resolvent G{f3) can be calcu- 
lated with linear scaling complexity for sufficiently large 



and sparse systems using several techniques, such as 
the Schulz iteration [^, the sparse approximate inverse 
[591 113 , and other methods • In a self-consistent cal- 
culation, where the Hamiltonian is changed in each iter- 
ation, or in a quantum molecular dynamics simulation, 
where both the overlap and the Hamiltonian is modified, 
we can efficiently update the new initialization 

Xi(new) — (Hncw ~ /3Snow) (36) 

from the previous iteration. If Xi(oid) and Xi(ncw) are 
sufficiently close the following scheme, based on Schulz's 
method (5& | , rapidly converges to the new normalization: 

Yo =-''^l(old)i 

^71+1 — 2Y^^ Y^ (Hncw /^ncwSncw)^^^ , (^'^) 
Xi(now) = Gnow(/3) = lim„^oo Y„. 

In this way, the cost can be reduced by using Schulz's 
method as an efficient iterative refinement technique. 



2. Non- Orthogonal Trace Correcting Purification 

A non-orthogonal second order trace correcting purifi- 
cation scheme is given by 

Xi = G(/3)-[H-(£,„i„-l)S]-\ (38) 



X„+i — F„(X„) — X„ -I- (T„(/ — X„S)X„, (39) 

where 

an = Sign[iVe - iVocc] (40) 



iVocc = Tr(SX„)] (41) 

and 

P = lim X„. (42) 

n — >oo 

Note that any of the initializations in Eqs. (|32|l - 134|l can 
be used to calculate Xi , realizing different levels of per- 
formance. 



B. Non-Orthogonal Perturbation Theory 

With an efficient normalization scheme in hand, given 
by Eq. H34|l . generalization of the density matrix pertur- 
bation theory to a non-orthogonal formulation follows, 
constituting the central result of this paper. At finite 
order, it provides the framework for computation of ba- 
sis set dependent response properties, including magnetic 
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response and geometric energy derivatives. The non- 
orthogonal extension is also useful for density matrix ex- 
trapolation in geometry optimization [b^ . 

The non-orthogonal perturbation theory below is de- 
veloped in the context of second order trace correcting 
purification. However, the formalism can be based also 
on other purification methods, such as grand canoni- 
cal purification ^^O, 26, 54], canonical purification [20l |. 
higher order trace correcting schemes and their hybrids 
|2n, 1^ , implicit purification at finite temperatures 
|5q|. or matrix sign function expansions [63Ll64j . 



1. Exact expansion (infinite order) 

Assume a perturbation of the Hamitonian and overlap 
matrix, 

H = H(o) +H', 

s = s(") + s'. 



(43) 



By analogy with Eq. I|17|l . with Xn and replaced 
by the non-orthogonal sequence X„ and purification poly- 
nomials F,i(X„), we have the non-orthogonal perturba- 
tions 



Ai=Fo(H(")+H')L-Fo(H(o)) 



S(0) 



(44) 



These perturbations generates a recursive purification se- 
quence that can be expanded to all orders in A, 



x„ = x„(")-i-ax„W-fa2x„(2) + 



(50) 



The initial expansion of Xi can be calculated using any 
of the normalizations given by Eqs. However, 
using the Green's function approach in Eq. (|34|l makes 
the expansion of Xi particularly simple. The terms are 



xl") 

X« 



XI- 



(2) ^ 
.(3) _ 



where 
and 



G(") 

-G(")T(i)G(o), 

„G(o)T(2)G(o) + G(o)T(i)G(o)T(i)G(o), 

_Q(0)rp(3)Q(0) + G(°)T(1)G(°)T(2)G("' - 



G(o) = (h(")-/3S(°))-i. 



(51) 
(52) 

(53) 



This initialization to various or der, X^™', in Eq. (|SU de- 
rives from the Dyson series 



A„+i =F„(X„-)-A„)|g-F„(X„)|g(o, . (45) 

Here F(X)|^ denotes the non-orthogonal purification 
with the metrics adapted to the overlap matrix A. With 
Fo(H) initiated by a non-orthogonal normalization in 
Eqs. (1221- 1231) and F„(X„) by Eq. ^ we have 



ri+l 



2 An - U„, iVocc < iVe, 



(46) 



where 



U„ = A„S(X„ + A„) + X„(SA„ + S'X„). (47) 

The occupation iVocc — Tr(S^"-'X„). Sine the perturba- 
tion theory is based on a perturbed projection, i.e. the 
difference between the purification of the perturbed and 
unperturbed sequence, the covariant (or contravariant) 
form of the difference A„ is preserved. This holds true 
also for the finite perturbation expansion described be- 
low. At convergence the non-orthogonal perturbed den- 
sity matrix is 

P = p(0) + lim A„. (48) 



2. Finite order perturbation 

Assume a perturbation expansion of the Hamiltonian 
and the overlap matrix, where 

H = H(0) + AH(i) + A2h(2) + . . . , 
S =SW+AS(1)+A2S(2) + ... . 



(49) 



G = G(°) - G(°) 



g(0) 



(54) 



from which terms in G can be collected order by order in 
A. A generalization to any order is straightforward. After 
the initialization of X^™^ we have for m — mmax, "T-max — 
1,...,1,0: 

r x(:^s(-')x«, N,cc>N, 

-v-("0 J i+j+k=m 

V i-\-j-\-k—7n 

(55) 

The sum is taken over all combinations (0, 1, ... , m) of 
i,j and k such that i + j + k = m. The occupation 
iVocc = Tr(S(")xl°)). At convergence the density matrix 
derivatives are given by 



1 9"P 



to! dX^ 



P(") = lim Xi™\ 



(56) 



A=0 



such that the density matrix perturbation expansion in 
a non-orthogonal representation is 



P-p(")+AP(i)+A2pC 



2) 



(57) 



This method for the calculation of density matrix re- 
sponse, including perturbations in the overlap matrix for 
a non-orthogonal representation, composes the central re- 
sult of this paper. 



IV. EXAMPLE 

To illustrate the non-orthogonal perturbation theory 
we have chosen a diatomic hydrogen io n H described in 
a basis set of two hydrogenic Is-orbitals |65l |. The overlap 
matrix S(i?) as a function of inter- atomic distance R (in 
units of Bohr radius oq) is given by 



51.1 — 82,2 — 1 

51. 2 = 82,1 = (1 



i?+(l/3)i?2 



(58) 



The matrix elements of the Hamiltonian H(i?) are 



Hl,l — H2,2 — Eq - 

Hi, 2 = H2,i = (Ef) 



i?-i(l-(l + i?)e-2fi) 



-kR-^ 
R)e-^, 
(59) 

where n — e^j (47reo) (set to 1 in the calculation). By ex- 
panding H and S in r = (i? — i?o) around the equilibrium 
distance i?o (or any other point) we have 



H = H(o) + rH(i) r2H(2) + 
S =S(o)+rSW+r2s(2)+.. 



where 



H(") = (TO!)-i5™H/ai?™ at i? = i?o, 
S(") = (TO!)-ia™S/ai?™ at R = i?o. 

(m) 



(60) 



(61) 



The initial perturbations A^'"^ are given by Eq. H51() 
and the recursive expansion is calculated as described in 
Eq. H55|l . At convergence the normalized density matrix 
derivatives P^*") are given by Eq. (|56|l . The expansion of 
the energy is given by collecting the energy 



£; = Tr 



(hW + rH(i) + . . .)(P(°) + rP(i) + ...)] (62) 



in orders of r = {R~ Ro)- Figure [21 shows the interaction 
potential E{R) as a function of inter-atomic distance in 
comparison to perturbation expansions up to 4**^ order 
at the equilibrium distance and at 4**^ order at a non- 
equilibrium inter-atomic distance. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have shown how density matrix per- 
turbation theory based on recursive purification can be 
generalized to include basis-set dependent perturbations. 
This makes it possible, for example, to calculate struc- 
tural response properties using local atomic-centered or- 
bitals within a reduced complexity formalism. Some key 
features of importance are: (1) an orbital-free density 
matrix formulation, which avoids the calculation of eigen- 
functions and eigenvalues, (2) very high order, mono- 
tonically increasing analytic approximation of the step 
function, (3) initial normalization of the Hamiltonian to 
fulfill the non-orthogonal commutation relation, which is 
preserved after each purification, and (4) the ability to 



' I ' I ' I ' 

— Exact Interaction Potential 
■ — Harmonic 2nd order 

— Anharmonic 3rd order 




Molecule 



Interatomic distance (R) 

FIG. 2: The analytic expansion of the energy E{R), Eq. (1621 . 
as a function of inter-atomic distance R for the molecule, 
using the non- orthogonal density matrix perturbation theory 
up to 4"' order. 



collect perturbations recursively, exactly (infinite order) 
or to any finite order, at each level of purification. 

A practical generalization of the Green's function ini- 
tialization in Eq. (|35|l is given by 



Xi = G(z) 



[(H-(en 



(63) 



which is stable for all z > 0. The value of z can be tuned 
to improve convergence and computational efficiency by 
optimizing the size of the band gap of the normalized 
spectra of Xi . The purification expansion is stable with 
respect to a complex generalization and the constant z 
can be extended to regions of the complex plane, in anal- 
ogy to Green's functions for complex energies. 

If an ill-conditioned non-orthogonal basis set is used we 
may run into numerical problems if we chose to trans- 
form the generalized eigenvalue problem to an orthog- 
onal representation. With the present formulation for 
non-orthogonal purification and perturbation theory, this 
congruence transform is avoided. Instead it is replaced by 
the calculation of G(z). However, if the condition num- 
ber of 6(2:) is smaller compared to the condition number 
of S, the numerical accuracy is improved. In addition, the 
back-transform from the orthogonal density matrix rep- 
resentation to the atomic orbital representation, which 
is necessary to calculate the electronic density expressed 
in the atomic orbital basis, is avoided within a purely 
non-orthogonal formalism. 

The example for the HJ molecule illustrates the ex- 
tension of the orbital-free density matrix perturbation 
theory to non-orthogonal representations. We have also 
applied the non-orthogonal method to recalculate the po- 
larizability of molecular clusters with results identical to 
previous calculations pol l57| . Since only matrix-matrix 
operations are used, the computational cost scales lin- 
early with system size for sufficiently large non-metallic 
systems, as was shown previously for the perturbed pro- 
jection scheme in an orthogonalized representation [29j. 
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The non-orthogonal density matrix perturbation theory 
can therefore efficiently be applied in calculations of re- 
sponse properties with a perturbation dependent basis 
set for large complex systems. 
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